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ABSTRACT 


This  study  compares  three  single  variable  search  methods  — 
Golden  Section,  cubic  interpolation  and  quadratic  interpolation. 
The  SUMT  nonlinear  program  was  used  for  the  comparison.   The 
OPT  subroutine  which  performs  the  single  variable  search  in 
SUMT  currently  uses  the  Golden  Section  method.   Two  different 
OPT  subroutines  were  written  which  implemented  cubic  inter- 
polation and  quadratic  interpolation.   Seven  test  problems 
which  contained  9-100  variables  and  2-20  constraints  were  used. 
The  comparison  was  made  on  computation  time  per  single  variable 
search  for  the  three  methods  and  the  number  of  function  evalua- 
tions per  single  variable  search  for  the  Golden  Section  and 
quadratic  interpolation  methods. 

A  single  variable  search  by  Lasdon,  Fox  and  Ratner  and 
one  by  Fletcher  and  McCann  were  also  discussed. 

The  results  showed  that  the  quadratic  interpolation  was 
slightly  faster  than  the  other  two  methods  and  required  fewer 
function  evaluations  per  single  variable  search  than  the 
Golden  Section  method.   Time  per  single  variable  search  was 
approximately  the  same  for  the  cubic  interpolation  and  Golden 
Section  methods.   Cubic  interpolation  required  fewer  points  to 
be  evaluated  than  the  other  two  methods  but  the  need  for 
gradient  evaluations  proved  to  be  costly  in  terms  of  computation 
time  per  single  variable  search. 
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I.   INTRODUCTION 

The  concept  of  optimization  has  grown  to  play  a  major 
role  in  the  analysis  of  many  complex  decision  and  allocation 
problems.   The  quality  of  the  analysis  not  only  depends  upon 
the  skill  and  good  judgement  used  in  interpreting  and 
modeling  the  problem  but  also  upon  the  reliability  and 
efficiency  of  the  vehicle  used  for  finding  a  solution  to  the 
problem.   The  solving  of  nonlinear  programming  problems  has 
seen  substantial  development  during  the  past  twenty  years 
and,  no  doubt,  will  increase  in  acceptance  and  popularity  in 
the  future . 

Within  the  Department  of  Defense  several  problems  of 
significant  importance,  such  as  weapons  targeting  and  alloca- 
tion, aircraft  and  power  plant  design  and  manpower  planning, 
to  name  a  few,  have  arisen  which  are  nonlinear  in  nature. 
Some  examples  of  industrial  nonlinear  programming  problems 
include  oil  refining,  curve  fitting,  inventory  and  logistics. 

Many  methods  exist  for  solving  nonlinear  programming 
problems.   Some  of  these  methods,  called  penalty  function 
methods,  are  based  on  transforming  a  constrained  minimization 
problem  into  a  sequence  of  unconstrained  minimization 
problems  whose  solutions  approach  the  solution  of  the 
original  constrained  problem. 

In  the  penalty  function  method,  as  well  as  in  most  other 
nonlinear  programming  algorithms,  a  one  dimensional 


minimization  search  is  used  repeatedly  in  the  solution 
process.   In  this  paper  several  different  one  dimensional 
minimization  search  methods  are  compared  in  the  context  of 
penalty  function  algorithms. 

The  different  one  dimensional  minimization  search  methods 
(also  called  single  variable  searches)  are  compared  by  the 
amount  of  computation  time  and  by  the  number  of  function, 
gradient  and  Hessian  evaluations  of  the  objective  function 
and  constraints  needed  to  find  a  solution  to  the  original 
constrained  optimization  problem. 

Section  II  of  this  paper  explains  the  formulation  of 
the  penalty  function  method  and  where  the  single  variable 
search  is  used  in  this  method.   The  nonlinear  computer 
program  SUMT  is  used  in  comparing • the  different  single 
variable  search  methods  and  is  also  explained. 

Section  III  formulates  the  different  single  variable 
search  methods  —  Golden  Section,  cubic  interpolation  and 
quadratic  interpolation  —  and  also  discusses  several  other 
methods  which  were  not  programmed  because  of  their 
complexity. 

Results  and  comparisons  of  the  different  methods  are 
stated  in  Section  IV  with  conclusions  and  suggestions  for 
further  study  in  Section  V.   The  Appendix  contains  the 
computer  codes  for  the  three  methods  that  were  programmed. 


II.   THE  PENALTY  FUNCTION  METHOD 

A.  THE  OPTIMIZATION  PROBLEM 

The  single  variable  search  is  a  major  part  of  finding 
an  optimal  solution  to  a  constrained  optimization  problem 
by  the  penalty  function  method.   Dr.  W.C.  My lander,  an 
author  of  the  SUMT  computer  code,  estimated  that  forty  percent 
of  the  computation  time  required  to  find  an  optimization 
problem  solution  using  SUMT  involved  single  variable  searches. 

The  constrained  optimization  problem  to  be  solved  is  of 
the  form: 

minimize:        f(x) 

subject  to:      g .  (x)  >_   0     j  =  l,...,m 

h.(x)  =  0     j  =  m+l,...,m+p  , 

where  x  is  an  n-dimensional  vector,  the  functions  f,  g  and  h 
are  continuous  and  may  be  either  linear  or  nonlinear,  and 
the  set  of  values  satisfying  the  inequality  constraints  has 
a  non-empty  interior,  i.e.  {x:  g. (x)  >  0,  j=l,...,m}  . 

B.  THE  SUMT  COMPUTER  PROGRAM 

The  computer  program  used  to  compare  the  different  single 
variable  search  methods  was  SUMT  -  Version  4  (Sequential 
Unconstrained  Minimization  Technique  for  Nonlinear  Programming) 
coded  in  FORTRAN  IV  by  W.C.  My lander,  R.L.  Holmes  and  G.P. 
McCormick  for  the  Research  Analysis  Corporation  [Ref.  1]. 
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It  was  written  to  experiment  with  different  methods  of 
solving  nonlinear  programming  problems. 

A  Guide  to  SUMT  -  Version  4  [Ref.  2]  is  the  program  manual 
which  contains  detailed  information  about  the  development  of 
the  code,  the  subroutines ,  options,  input-output,  limitations, 
user-supplied  subroutines,  data  deck  setup  and  a  detailed 
example  of  the  use  of  the  program  to  solve  a  nonlinear 
programming  problem. 

The  method  used  by  SUMT  to  solve  the  original  constrained 
optimization  problem  is  described  in  detail  in  Chapter  8  of 
the  book  on  nonlinear  programming  by  Fiacco  and  McCormick 
[Ref.  3] .   Basically,  SUMT  solves  the  constrained  problem 
by  finding  the  solutions  of  a  sequence  of  unconstrained 
problems  which  approach  the  solution  of  the  original  con- 
strained problem. 

Previous  versions  of  SUMT  used  different  penalty  functions 
for  approximating  the  solution  to  the  constrained  problem  but 
the  function  currently  used  to  sequentially  solve  the  problem 
is  of  the  form: 


m  m+p         2 

P(x,r)  =  f(x)  -  r   Z   In  g.  (x)   +    E    [h.  (x)JVr 

j=l     J  j=m+l    J 


where  the  parameter  r  determines  the  severity  of  the  penalty 
and  consequently  how  closely  the  unconstrained  problem 
solution  approaches  the  solution  to  the  constrained  problem. 

The  penalty  function  of  SUMT  uses  the  mixed  interior 
point-exterior  point  method  described  in  detail  in  Chapter  4 
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m 

of  Ref.  2.   The  term,   r  Z   In  g.(x)  ,  involving  the 

j-1  D 

inequality  constraints  is  from  the  interior  point  method 

which  is  also  known  as  the  barrier  method.   As  the  x  vector 

approaches  the  boundary  of  the  infeasible  region  for  these 

constraints,  this  term  approaches  infinity.   The  term, 
m+p         2 
E    [h . (x) ]  /r  ,  involving  the  equality  constraints  is 
j=m+l   J 
from  the  exterior  point  method.   As  the  x  vector  moves 

farther  away  from  the  feasible  region  for  these  constraints, 

this  term  approaches  infinity.   The  behavior  of  these  two 

terms  is  shown  graphically  in  Figure  1.   How  much  of  a 

penalty  these  terms  add  to  the  penalty  function  depends  upon 

the  value  of  the  parameter  r. 

The  solution  procedure  requires  minimization  of  P(x,r) 
over  those  x  satisfying  the  conditions   g.(x)  >  0  ,  j=l,...,m 
for  r  =  r,,r2/...  where  r,  <  r~  <...  <  r,  <  ...  <  0  .   Under 
mild  conditions  on  the  original  NLP  problem  the  minima  of 
the  sequential  unconstrained  problems  approach  the  solution 
of  the  original  constrained  problem  as  r.  +  0.   The  condi- 
tions to  be  met  for  the  method  to  solve  the  original  problem 
are  described  in  detail  in  the  SUMT  program  manual  [Ref.  2]. 

When  the  NLP  problem  is  convex  the  SUMT  program  also 
generates  a  sequence  of  points  that  are  feasible  for  an 
optimization  problem  called  the  dual  of  the  original  prob- 
lem.  The  maximum  function  value  of  the  dual  feasible  points 
and  minimum  function  value  of  the  unconstrained  problem 
bracket  the  optimal  solution  of  the  constrained  problem. 
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P(x,r) 


S      INFEASIBLE 


a.   INTERIOR  POINT  METHOD  INVOLVING  INEQUALITY  CONSTRAINTS 


P(x,r) 


b.      EXTERIOR  POINT  METHOD  INVOLVING  EQUALITY  CONSTRAINTS 


FIGURE    1.       BEHAVIOR   OF    THE    PENALTY    TERMS 
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As  the  unconstrained  problem  is  minimized  and  the  number  of 
dual  feasible  points  found  increases  the  difference  between 
the  two  solutions  decreases  and  the  optimal  solution  is  more 
accurately  approximated. 

SUMT  makes  use  of  the  theory  derived  in  Fiacco  and 
McCormick's  book  [Ref.  3,  Chapter  5]  that  uses  the  Lagrange 
extrapolation  technique  to  approximate  the  solution  of  the 
constrained  problem  when  more  than  one  minimum  of  the 
unconstrained  problem  has  been  found.   By  using  these 
extrapolation  approximations  the  solution  converges  faster 
and  achieves  a  closer  approximation  to  the  actual  solution 
than  the  minima  of  the  unconstrained  problems. 

The  logic  of  the  computer  program  taken  from  the  SUMT 
manual  [Ref.  2]  follows. 

STEP  1. 

Find  an  initial  starting  point  x  within  the  interior 
feasible  region  such  that  g . (x)  >  0  ,  j=l,...,m  .   If  such 
an  initial  point  is  not  given  or  is  not  feasible,  the  program 
uses  the  optimization  method  to  find  one. 

STEP  2. 

Determine  r, ,  the  initial  value  of  r,  which  can  either 
be  an  input  parameter  or  a  rule  for  choosing  r1  can  be 
specified. 

STEP  3. 

Determine  the  minimum  of  the  unconstrained  penalty 
function  for  the  current  value  of  r. 
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STEP  4. 

Estimate  a  better  solution  using  the  extrapolation 
technique,  if  possible. 

STEP  5. 

Computation  is  terminated  if  the  stopping  criteria  are 
satisfied  thus  yielding  an  approximate  solution  to  the 
original  problem  with  accuracy  depending  upon  the  stringency 
of  the  stopping  criteria. 

STEP  6. 

If  the  stopping  criteria  are  not  satisfied  select 

rk+l  K   rk  * 
STEP  7. 

Go  to  Step  3. 

A  simplified  flow  diagram  of  the  computer  logic  taken 
from  the  SUMT  manual  [Ref.  2,  p.  5]  is  shown  in  Figure  2. 

The  single  variable  search  problem  is  involved  in  Step  3 
of  the  procedure.   In  order  to  minimize  the  unconstrained 
penalty  function  the  program  chooses  a  search  direction  in 
n-space.   It  is  believed  that  by  searching  in  this  direction 
from  the  starting  point  of  the  subproblem  the  penalty  function 
will  initially  decrease.   A  one  dimensional  search  is  then 
made  along  that  direction  until  a  local  minimum  is  found. 
At  the  solution  to  this  one  dimensional  search  problem  a 
new  search  direction  is  computed  and  a  new  one  dimensional 
search  is  performed.   The  one  dimensional  searches  are 
repeated,  each  time  with  a  new  search  direction,  until  the 
unconstrained  problem  for  a  given  r  is  solved. 
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FIGURE    2.       SUMT   Program   Flow   Diagram    [Ref.    2,    p.    5] 
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SUMT  allows  the  user  to  choose  one  of  three  procedures 
to  determine  the  search  direction  —  by  Newton's  Method  using 
first  and  second  partial  derivatives  of  the  unconstrained 
penalty  function,  by  using  the  negative  of  the  gradient  of 
the  penalty  function  or  by  a  variable  metric  method  which  is 
taken  from  an  algorithm  of  the  Fiacco  and  McCormick  book 
[Ref.  3,  pp.  170-175]. 

Finding  a  solution  to  the  unconstrained  n  dimensional 
problem  thus  requires  several  single  variable  searches. 
Therefore,  the  more  efficient  the  single  variable  search  is, 
the  faster  the  computation  time  will  be.   Given  a  direction 
of  search,  the  SUMT  program  uses  the  subroutine  OPT  to 
perform  the  single  variable  search.   In  this  paper,  the 
presently  used  single  variable  search  method  in  the  OPT 
subroutine  is  compared  to  other  single  variable  search  methods 
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III.   SINGLE  VARIABLE  SEARCH  METHODS 

Since  a  major  portion  of  finding  a  solution  to  the 
constrained  optimization  problem  involves  single  variable 
searches,  the  SUMT  program  was  used  to  compare  several 
single  variable  search  methods. 

The  current  OPT  subroutine  of  SUMT  uses  the  Golden 
Section  search  method  to  perform  the  single  variable  search. 
Two  different  single  variable  search  methods  were  programmed 
to  replace  the  current  OPT  subroutine.   The  first  method 

uses  cubic  interpolation  once  the  local  minimum  is 

bracketed  to  accelerate  convergence  to  that  local  minimum 

and  the  second  method  programmed  uses  quadratic  interpolation 

to  accelerate  convergence  once  three  feasible  points  that 

bracket  the  local  minimum  are  found.   Single  variable  search 

methods  by  Lasdon,  Fox  and  Ratner  [Ref.  4]  and  Fletcher  and 

McCann  [Ref.  5]  are  also  discussed. 

The  object  of  the  single  variable  search  is  to  find  a 

scalar  9  called  the  step  length  such  that  x*  =  x  +  9s 

o 

where  x*  is  the  x  vector  that  corresponds  to  the  local 
minimum  of  the  penalty  function  along  the  search  direction  s 
from  an  initial  starting  point  x  .   For  ease  of  notation 
penalty  function  values  are  shown  as  a  function  of  step 
length  9  in  each  single  variable  search  since  the  initial 
starting  point  x  ,  the  search  direction  s  and  the  parameter 
r  are  constant,  i.e.  P(9)  =  P(x  +9*s,r). 
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The  optimal  step  length  could  possibly  be  found  by 
moving  along  the  search  direction  in  a  random  manner. 
However,  this  would  not  be  a  very  logical  way  to  attack  the 
problem  and  would  lead  to  many  needless  evaluations  of  the 
penalty  function  which  would  increase  computation  time. 

The  OPT  subroutine  is  supplied  with  a  search  direction 
and  an  initial  point  from  which  to  search.   It  must  call 
other  subroutines  of  the  SUMT  program  for  functional, 
gradient  or  Hessian  evaluations  and  is  expected  to  return 
with  a  penalty  function  value  and  corresponding  x  vector 
which  is  the  local  minimum  for  that  single  variable  search. 
It  is  assumed  that  the  penalty  function  has  a  finite  local 
minimum  along  the  given  search  direction. 

The  following  subsections  formulate  and  explain  the 
different  single  variable  search  methods. 

A.   GOLDEN  SECTION  SEARCH  METHOD 

The  Golden  Section  search  method  uses  the  limiting 
properties  of  the  Fibonacci  series  [Ref .  6] .   Fiacco  and 
McCormick's  book  [Ref.  3]  define  the  steps  of  the  procedure 
used  in  the  Golden  Section  method. 

STEP  1. 

First  an  upper  bound  step  length  0__  is  obtained  with  the 
lower  bound  step  length  8   initially  equal  to  zero.   9y  is 
obtained  by  evaluating  the  penalty  function  at  successive 
step  lengths  which  are  in  the  limiting  Fibonacci  ratio, 
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k 
(1  +  /D/2    z    1.618    ,      that   is,    Q^  =      E       (1.618)1    ,      where 

i=0 

k  is  the  smallest  integer  j  such  that 

P  [  I       (1.618) *]  >  P  [  I       (1.618)*] 
£=0  £=0 

j"2         £ 
with  9_  updated  as   E   (1.618)   . 

L  £=0 

STEP  2. 

The  interval  bounding  9*,  the  optimal  step  length ,    is 
reduced  by  computing  two  other  step  lengths  within  the 
interval  (6_,9„).   Their  corresponding  values  are: 


)L  +  0.382(9U  -  9L) 


>L  +  0.618(9  -  9  ) 


STEP  3. 

The  penalty  function  values  of  the  two  interior  step 
lengths,  9   and  9,  ,  are  then  compared. 

STEP  4. 

If  P(9  )  <  P(9K),  then  the  optimal  step  length  which 
a       o 

corresponds  to  the  local  minimum  of  the  single  variable 

search  should  be  in  the  interval  (8L,6,  )  if  the  penalty 

function  is  unimodal.   Because  of  the  fact  that 

0.382/0.618  *  0.618,  reassign  Q^   =  9fa,  9fa  =  9a  and  calculate 

a  new  9^  as  computed  in  Step  2.   Return  to  Step  3. 
a 
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STEP  5. 

If  P(8a)  >  p(9b)/  then  the  optimal  step  size  should  be 

in  the  interval  (9  ,9TT).   Reassign  6T  =  9  ,  9   =9,  and 

a   u  Jj    a   a    d 

calculate  a  new  9fc  as  computed  in  Step  2.   Return  to  Step  3. 

STEP  6. 

If  the  penalty  function  values  at  both  interior  step 
lengths  are  equal,  reassign  9L  =  9  ,  Q^   =  9fa  and  calculate 
a  new  9   and  9^  as  computed  before  by  returning  to  Step  2. 

STEP  7. 

When  the  stopping  criteria  is  satisfied  as  explained 
in  subsection  III.D,  9*  is  approximated  by  (9L  +  9u)/2 
yielding  an  x  vector  x*  =  x  +  9*S  which  corresponds  to  an 
approximated  local  minimum  P(9*)  of  the  single  variable 
search. 

A  flow  diagram  of  the  Golden  Section  method  presently 

used  in  the  OPT  subroutine  is  shown  in  Figure  3.   Certain 

parts  of  the  computer  code  that  were  the  same  in  each  of 

the  OPT  subroutines  are  discussed  in  subsection  III.D. 

The  Golden  Section  OPT  subroutine  initially  updates  9 

j_i 

as  it  finds  more  feasible  points  by  increasing  the  step 
length  until  it  finds  a  point  which  is  infeasible  or  where 
the  penalty  function  value  is  larger  than  the  penalty 
function  value  at  the  previous  feasible  point.   If  it  finds 

an  infeasible  step  length,  the  penalty  function  is  assigned 

3  6 
a  very  high  value  (10   ) ,  the  next  to  the  last  feasible 

step  length  is  assigned  9-  and  the  last  feasible  step  length 

assigned  as  the  lower  interior  step  length  9  .   If  the  first 
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FIGURE  3.   Golden  Section  OPT  Subroutine 
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step  length  from  the  initial  point  is  infeasible,  the  lower 
interior  step  length  9a  =  0.382  9  ,  since  only  one  feasible 
step  length  (9T  =  0)  has  been  found.   OPT  then  uses  8„  and 

1j  U 

8   to  compute  the  upper  interior  step  length 
9b  =  9a  +  0.382 (9n  -  9a) .   The  optimal  step  length  is 
bracketed  and  the  subroutine  continues  to  reduce  the 
interval  of  uncertainty  by  comparing  penalty  function  values 
of  the  interior  step  lengths  until  it  finds  two  values  or 
an  interval  that  satisfies  the  stopping  criteria. 

In  finding  the  minimum  of  each  single  variable  search, 
the  Golden  Section  method  only  requires  penalty  function 
evaluations  which  are  compared  to  show  where  the  penalty 
function  decreases  to  a  local  minimum  and  then  increases  to 
infinity  as  the  infeasible  region  boundary  is  approached. 
The  factor  by  which  the  interval  of  uncertainty  is  reduced 
remains  constant  in  this  method.   Therefore,  any  information 
about  the  behavior  of  the  penalty  function  other  than  the 
fact  that  the  optimal  step  length  is  somewhere  within  the 
interval  of  uncertainty  is  disregarded. 

B.   CUBIC  INTERPOLATION  METHOD 

The  cubic  interpolation  method  was  programmed  for 
comparison  with  the  Golden  Section  and  quadratic  interpola- 
tion methods.   Basically,  this  method  approximates  the 
penalty  function  by  a  cubic  fit  and  then  finds  the  minimum 
of  the  cubic  equation  to  approximate  the  optimal  step 
length  of  the  single  variable  search. 
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In  order  to  approximate  the  penalty  function  by  a  cubic 
fit  in  the  single  variable  search,  two  feasible  step  lengths 
which  bracket  the  local  minimum  along  with  their  directional 
derivative  values  are  needed.   This  method  actually  consists 
of  two  distinct  parts  —  the  bracketing  section  to  find  the 
two  feasible  points  and  the  cubic  interpolation  section  to 
find  a  minimum  of  the  cubic  function.   The  single  variable 
search  method  is  supplied  with  an  initial  starting  point  x  , 
the  penalty  function  value  at  x   [P(0)],  the  gradient  of 
the  penalty  function  at  x   [VP(0)]  and  a  search  direction  s. 

The  steps  of  the  cubic  interpolation  method  follow. 

STEP  1. 

Since  the  directional  derivative  at  the  starting  point 
P1  CO)  is  negative,  by  finding  a  step  length  that  has  a 
corresponding  positive  directional  derivative,  the  local 
minimum  along  the  given  direction  of  search  will  be  bracketed 
if  the  penalty  function  is  unimodal. 

The  step  length  is  increased  to  find  an  upper  step 
length  9   until  an  n  is  found  such  that 

n 
P'  [(  Z   21)  -  1]  >  0   . 
i=l 

If,  while  initially  increasing  the  step  length,  an  infeasible 
point  is  encountered,  the  increment  by  which  the  last  step 
length  was  increased  is  halved  and  subtracted  from  the 
current  step  length.   This  is  continued  until  a  step  length 
that  is  feasible  is  found.   If  its  directional  derivative  is 
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negative,  the  increment  is  increased  by  one  half  the 
increment  length  and  added  to  the  current  step  length. 
This  is  continued  until  a  step  length  with  a  positive 
directional  derivative  is  found. 

The  lower  step  length  0,  is  always  the  last  feasible 
step  length  encountered  with  a  negative  directional 
derivative. 

STEP  2. 

Once  two  step  lengths  have  been  found  that  bracket  the 
local  minimum  along  the  search  direction,  the  cubic  equation 
is  used  to  compute  a  step  length  which  corresponds  to  the 
minimum  of  the  cubic  fit.   This  step  length  is  found  by  the 
equation: 


p'<eD)  +  u2  -  u1 
e  -  eTT  -  (eTT  -  eT )    [  ^ 


u  1/         P'  (9    )    -   P*  (8L)    +    2U. 


where  .  p(eT  -  p(9tt) 

ox-  P'(eL)   +  p'(9u)   -  3   [     e^  .  6u      1 


u2  =   [u12  -  p'te^p'Ojj)  ]1/2 


STEP    3. 

P'(9)  is  evaluated  and  if  it  is  negative  then  reassign 
=  9  and  the  left  side  is  discarded.   If  P*(8)  is  positive 


L 


then  reassign  9n  =9   and  the  right  side  is  discarded. 
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STEP  4. 

If  the  stopping  criteria  is  satisfied,  the  two  penalty 
function  values,  P(8)  and  P(9  )  or  P(8)  and  P(6  ),  depending 
upon  P'(0),  are  compared  and  the  smallest  one  is  returned 
as  the  local  minimum  P(0*)  of  the  single  variable  search 
with  its  corresponding  x  vector,  x  +9*s. 

STEP  5. 

If  the  stopping  criteria  is  not  satisfied  return  to 
Step  2. 

A  flow  diagram  of  the  cubic  interpolation  OPT  subroutine 
is  shown  in  Figure  4.   The  parts  of  this  subroutine  that  are 
the  same  in  the  other  OPT  subroutines  are  discussed  in 
subsection  III.D. 

The  cubic  interpolation  method  uses  more  information  at 
each  feasible  step  length  to  find  the  local  minimum  of  the 
single  variable  search  than  the  Golden  Section  or  quadratic 
interpolation  methods.   It  requires  a  gradient  evaluation 
for  the  penalty  function  at  each  feasible  step  length  in 
order  to  compute  the  directional  derivative.   A  gradient 
evaluation  may  require  many  function  evaluations  depending 
upon  the  number  of  variables  and  constraints  in  the  original 
constrained  problem.   However,  with  this  extra  information 
the  local  minimum  should  be  found  by  evaluating  fewer  step 
lengths  in  the  cubic  interpolation  section  of  this  method. 
As  the  value  of  the  parameter  r  decreases  with  each 
unconstrained  problem,  the  penalty  function  rises  more 
abruptly  as  the  infeasible  region  boundary  is  approached. 
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FIGURE    4.      Cubic    Interpolation   OPT    Subroutine 
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This  makes  it  more  difficult  to  find  a  step  length  with  a 
positive  directional  derivative  which  is  required  to  bracket 
the  minimum. 

C.   QUADRATIC  INTERPOLATION  METHOD 

A  quadratic  interpolation  single  variable  search  was  also 
programmed  for  comparison  with  the  two  previously  discussed 
methods . 

This  method  requires  three  feasible  step  lengths  and 

their  penalty  function  values  which  obey  the  relationship 

P(0_)  >  P(9J  <  P(8TT)  where  9_  <  9   <  9TT  .   With  this 
L        M        U  L     M     U 

information  a  quadratic  fit  is  made  with  its  minimum 
approximating  the  minimum  of  the  penalty  function  for  the 
given  search  direction.   The  penalty  function  is  evaluated 
at  the  interpolated  step  length  which  minimizes  the  quadratic 
function  and  compared  to  P ( 9  )  allowing  the  interval  of 
uncertainty  to  be  decreased.   The  process  is  repeated  until 
the  minimum  of  the  quadratic  fit  approximates  the  local 
minimum  of  the  penalty  function  or  the  interval  of  uncertainty 
becomes  small  enough  to  satisfy  the  stopping  criteria. 

This  method  also  consists  of  two  sections  —  the  bracketing 
section  and  the  quadratic  interpolation  section.   The  single 
variable  search  is  supplied  with  an  initial  x  vector  x  ,  its 
penalty  function  value  P(0)  and  a  search  direction  s. 

The  steps  of  the  quadratic  interpolation  search  method 
follow. 
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STEP  1. 

The  step  length  is  increased  from  the  initial  starting 
point  to  find  an  upper  step  length  9   until  an  n  is  found 
such  that 


P[  (  Z   21)  -  1]  >  P[  (  Z  21)    -  1]   . 
i=l  i=l 


If,  while  initially  increasing  the  step  length,  an  infeasible 
step  length  is  encountered,  the  increment  by  which  the  last 
step  length  was  increased  is  halved  and  subtracted  from  the 
current  step  length.   This  is  continued  until  a  feasible  step 
length  is  found.   If  its  functional  value  is  less  than  P(9T), 

Li 

the  increment  is  increased  by  one  half  the  increment  length. 
This  is  continued  until  a  feasible  step  length  is  found  that 
has  a  larger  penalty  function  value  than  the  previous  step 

and 
9M  =  9  ,  always  ensuring  that  POt)  >  P ( 9«) • 

If  only  one  feasible  step  length  is  found  in  obtaining 
the  upper  step  length,  where  in  this  case  P(9y)  >   p(eL) 
(which  indicates  the  minimum  is  bracketed) ,  the  last 
increment  is  halved  and  subtracted  from  the  upper  step 
length  which  yields  the  interior  step  length  6m. 

STEP  2. 


length.   Before  this  is  repeated,  reassign  9L  =  9M 


Now  that  the  local  minimum  is  bracketed,  the  quadratic 
function  is  used  to  find  a  step  length  where  the  derivative 
of  the  quadratic  fit  vanishes. 
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_  1   b23P(9L)  +  b31P(V  +  b12P(V 
2   a23P(9L)  +  a31P(6M)  +  a12P(9u) 

where  a..  -9.-0.  ,   b..  =  8.2  -  9.2  , 

i,j  =  1,2,3  where  1  =  L,  2  =  M,  and  3  =  U. 

STEP  3. 

If  9  >  6M,  then  the  left  side  is  discarded  and  reassign 

9L  =  9M  and  9M  =  9*   If  9  <  9M'  then  the  ri9ht  side  is 

discarded  and  reassign  9TT  =  9..  and  9„  =  9. 

U    M      M 

STEP  4. 

P(9M)  is  evaluated  and  if  the  stopping  criteria. are 
satisfied,  the  smallest  of  P(9L),  P(9M)  and  P(9  )  is  returned 
as  P(9*)  with  the  corresponding  x  vector,  x  +9*s. 

STEP  5. 

If  the  stopping  criteria  are  not  met,  return  to  Step  2. 

A  flow  diagram  of  the  quadratic  interpolation  method  OPT 
subroutine  is  shown  in  Figure  5. 

The  quadratic  interpolation  method  only  requires  function 
evaluations  to  perform  the  single  variable  search.   It  uses 
the  knowledge  of  three  feasible  step  lengths  and  their 
penalty  function  values  to  approximate  where  the  local 
minimum  of  the  single  variable  search  is  located.   However, 
it  must  find  an  upper  feasible  step  length  that  is  within  the 
interval  where  the  penalty  function  is  increasing  to  infinity 
and  also  this  feasible  step  length  must  have  a  penalty 
function  value  that  is  greater  than  that  of  the  interior 
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feasible  step  length.   This  causes  the  interval  containing 
possible  feasible  upper  step  lengths  to  be  reduced  even  more. 
As  the  parameter  r  decreases  for  each  unconstrained  sub- 
problem,  the  penalty  function  increases  more  abruptly  from 
the  minimum  as  the  infeasible  region  boundary  is  approached. 
This  also  increases  the  difficulty  in  bracketing  the  minimum. 

Once  the  three  feasible  step  lengths  are  found,  the 
quadratic  interpolation  can  converge  to  the  local  minimum 
thereby  solving  the  single  variable  search. 

D.   SECTIONS  COMMON  TO  THE  DIFFERENT  METHODS 

Some  parts  of  the  three  different  OPT  subroutines  were 

made  similar  so  that  they  could  be  compared  under  the  same 

conditions. 

Essentially  the  same  stopping  criteria  was  used  in  each 

subroutine.   When  the  ratios  of  the  two  x  vectors  that 

bracketed  the  minimum  or  ratios  of  their  penalty  function 

-7 
values  were  within  10    of  one,  the  subroutine  returned  to 

the  program  with  a  local  minimum  and  a  corresponding  x  vector 
as  the  solution  to  the  single  variable  search.   These  criteria 
were  checked  in  the  Golden  Section  method  whenever  two 
interior  step  lengths  had  been  computed  and  in  the  inter- 
polation methods  both  after  the  minimum  had  been  bracketed 
and  also  after  each  interpolated  step  length  had  been 
computed. 

In  the  three  different  methods,  the  subroutines  also 
contained  counters  which  would  cause  the  search  to  stop 
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after  ten  interpolations  in  the  cubic  interpolation  method, 
twenty  interpolations  in  the  quadratic  interpolation  method 
or  twenty-five  feasible  points  in  the  Golden  Section  method. 
In  the  test  problems  that  were  run  these  counters  were  never 
reached,  thereby,  never  satisfying  this  stopping  criterion. 

In  the  cubic  interpolation  method  another  stopping 
criterion  was  used  along  with  the  ratio  tests.   When  the 

cosine  of  the  angle  between  the  gradient  of  the  penalty 

-7 
function  and  the  search  direction  vector  was  less  than  10 

the  subroutine  returned  with  a  local  minimum  and  x  vector. 

It  was  stated  earlier  that  it  was  hoped  that  the  penalty 

function  would  initially  decrease  in  the  given  search 

direction  from  the  initial  starting  point.   In  the  three 

different  methods,  if  the  ratio  of  any  element  of  the 

current  x  vector  and  the  corresponding  element  in  the  initial 

-7 
x  vector  came  within  10   of  one,  the  single  variable  search 

was  restarted  in  the  opposite  search  direction.   In  the 

interpolation  methods  if  after  one  hundred  tries  the  minimum 

was  not  bracketed  or  the  step  length  had  been  decreased  more 

than  twenty  consecutive  times,  the  search  direction  was  also 

reversed  and  the  single  variable  search  was  restarted. 

The  three  different  OPT  subroutines  called  the  subroutine 

EVALU  for  penalty  function  evaluations.   EVALU  would  return 

back  a  variable  after  each  call  up  which  would  show  if  the 

step  length  was  feasible.   The  cubic  interpolation  method 

called  the  subroutine  GRAD  for  the  gradient  of  the  penalty 

function  when  it  was  needed  to  determine  the  directional 

derivative. 

33 


E.   OTHER  SINGLE  VARIABLE  SEARCH  METHODS 

Two  other  single  variable  search  methods  were  also 
looked  at  but  were  not  programmed  because  of  their  complexity 
and  their  need  for  changing  other  portions  of  the  SUMT 
program  than  just  the  OPT  subroutine. 
1.   Lasdon,  Fox,  and  Ratner  Method 

The  single  variable  search  developed  by  Lasdon,  Fox 
and  Ratner  [Ref.  4]  uses  a  penalty  function  that  doesn't 
contain  equality  constraints  and  the  penalty  term  for  the 
inequalities  is  also  different  than  the  SUMT  penalty  function 
The  unconstrained  minimization  problem  is: 


m     1 
minimize  P(xfr)  =  F(x)  +  r   E   - — 7 — r 

i=l   Gi(x) 


where  F  is  the  objective  function  and  the  G.  *s  are  the 
inequality  constraints.   This  model  could  be  modified  to 
handle  equality  constraints  and  a  penalty  function  like  that 
used  in  SUMT. 

The  single  variable  search  method  is  in  three  distinct 
stages  —  linear  approximation,  quadratic  approximation  and 
cubic  interpolation.   Linear  approximations  are  made  for  the 
objective  function  and  each  constraint  using  F(0),  F'(0), 
G.(0)  and  G.'(O).   A  step  length  (9-,)  which  minimizes  the 
penalty  function  consisting  of  the  linear  approximations  is 
then  calculated.   The  information  used  to  make  the  linear 
approximations  is  then  used  along  with  F(91),  F'O^,  G.  (01) 
and  G-'(9-|)  to  make  quadratic  approximations  of  the 
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objective  function  and  constraints.   The  penalty  function 
consisting  of  quadratic  approximations  is  then  minimized 
by  computing  92«   If  the  directional  derivative  of  the 
original  penalty  function  is  positive  at  92,  cubic  inter- 
polation is  applied  yielding  0*,  the  solution  to  the  single 
variable  search.   If  the  stopping  criterion  is  satisfied 
during  any  stage,  the  single  variable  search  is  terminated 
yielding  a  local  minimum  and  a  solution  to  the  search. 

The  following  steps  outline  the  Lasdon,  Fox  and 
Ratner  procedure. 

STAGE  1. 

The  objective  function  and  inequality  constraints 
of  the  original  constrained  problem  are  approximated  by 
using  the  given  values  at  the  initial  starting  point. 


f1(9)  =  F(0)  +  F1 (0) 


gj_.Ce)  =  G.(0)  +  G.'(0)      j  =  l,...,m. 


The  approximating  function  for  P(9)  is 


m 


1 


Pl(9)  =  fL(9)  +  r  .^^-TeT 

The  smallest  positive  zero  of  the  linear  approximations  of 
the  inequality  constraints  is  found  by  taking  the  smallest 
-G. (0)/G. ' (0)  over  all  j.   This  value  is  designated  the 
upper  step  length  9y. 

The  step  length  corresponding  to  the  minimum  of  the 
approximated  penalty  function  is  found  by  doing  four  or  five 
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Newton's  method  iterations  using  ■=■  6   as  the  starting  point 


1  ~  2  8U 1 

Pi" <?  9U> 


The  result  of  this  stage  is  a  step  length  that 
approximates  the  optimal  step  length  of  the  single  variable 
search.   If  the  stopping  criterion,  explained  following 
Stage  3,  is  not  satisfied  proceed  to  Stage  2. 

STAGE  2. 

The  same  procedure  is  performed  in  this  stage  as  in 
Stage  1  with  the  exception  that  quadratic  approximations  to 
the  objective  function  and  constraints  are  made  instead  of 
linear  approximations.   The  second  stage  is  entered  with 
values  of  F(0,)  and  G.  (0, )  along  with  the  information  used 
to  make  the  linear  approximations.   The  quadratic  approxima- 
tions to  the  objective  function  and  constraints  are 

f2(6)  =  ae2  +  be  +  c 

F(e1)  -  be1  -  c 

where      a  =  

b  =  F' (0)  ,   c  =  F(0) 


and 


g9.(0)  =  d.02  +  e.   +  f .     j  =  l,...,m 

G  (ex)  -  ej01  -  f. 
where      d  .  =  — * ~ 

3  91 

ej  =  Gj'tO)  ,   fj  =  G^O)  . 
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If  Gj(81)  is  infeasible  for  some  j,  then  the  smallest 
positive  root  over  all  j  of 


-e.  ±  /e.2  -  4d.f . 

—3 3 J  J    =  e 


2dj 

is  used  as  the  starting  step  length  for  using  Newton's 
method  to  minimize  the  approximated  penalty  function 

m      , 

p2(6)  =  f,(9)  +  r   S   W   ' 

2        2         .=1   g2j(9) 

The  step  length  which  minimizes  P?(e)  should  be 
found  after  four  or  five  iterations  of  Newton's  method. 

P  '(6  ) 

9=9-  -2 L. 

P2"(91) 

If  91  is  infeasible,  use  -j  9    in  Newton's  method. 
If  the  stopping  criterion  is  not  satisfied  after  finding 
92  continue  to  Stage  3. 

STAGE  3. 

Direct  cubic  interpolation  of  the  penalty  function 
is  used  in  this  stage.   If  P'(92)  is  positive,  the  minimum 
is  bracketed  and  cubic  interpolation  can  be  performed  to 
find  a  step  length  which  corresponds  to  the  minimum  of  the 
penalty  function.   If  the  minimum  has  not  been  bracketed, 
that  is,  if  P'(92)  <  0,  a  bracketing  procedure  like  the  one 
used  in  the  cubic  interpolation  method  of  subsection  III.B 
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must  be  performed  to  find  a  step  length  that  has  a  positive 
directional  derivative.   Once  this  step  length  is  found, 
cubic  interpolation  is  used  to  find  a  step  length  which 
minimizes  the  cubic  function.   If  the  new  step  length  does 
not  satisfy  the  stopping  criteria,  additional  cubic  inter- 
polations are  made  using  the  two  points  which  most  closely 
bracket  the  minimum  of  the  penalty  function. 

When  two  penalty  function  values  are  nearly  equal 
or  when  the  interval  between  step  lengths  that  bracket  the 
minimum  becomes  very  small,  an  optimal  step  length  is  deter- 
mined by 

Q*  =  P'(ea)eb  -  P'(eb)9a 

P'  (9  )  -  P1  (6.  ) 
a        b 

where  9   is  the  lower  step  length  and  0.  is  the  upper  step 
a  d 

length  (8   <  8*  <  9b) . 

The  stopping  criterion  used  in  this  method  is  the 

I sTVP I 
cosine  test   |  cos  (j>  |  =  ' ' ,      where  <£  is  the  angle 

I  I  s  |  |  |  | VP  |  | 
between  the  search  direction  vector  and  gradient  of  the 
penalty  function.   When  the  cosine  of  the  angle  is  less  than 
10   ,  the  stopping  criterion  is  satisfied. 

The  authors  found  that  in  using  this  method  on  test 
problems,  the  stopping  criterion  was  often  satisfied  at  the 
end  of  Stage  2.   In  their  comparisons  of  this  method  with  a 
cubic  interpolation  method,  they  found  that  it  performed 
the  single  variable  search  much  more  efficiently  [Ref.  4, 
p.  295]. 
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This  method  requires  gradients  of  the  objective 
function  and  constraints.   In  the  SUMT  program  these  values 
are  not  stored  so  implementing  this  method  would  require 
more  storage  or  more  gradient  evaluations  along  with 
changes  in  other  subroutines  and/or  addition  of  new 
subroutines . 

2 .   Fletcher-McCann  Method 

The  Fletcher-McCann  method  [Ref.  5,  pp.  210-212] 
uses  an  approximation  of  the  original  unconstrained  penalty 
function  to  find  the  local  minimum  of  the  single  variable 
search.   The  penalty  function  approximation  used  is  of  the 
form 

T(8)  =  a  +  be  +  c92  +  5  d  5  , 


U 

where  the  parameters  a,  b,  c,  d  and  0   are  determined  by 
using  objective  function  and  constraint  information  at  two 
feasible  step  lengths.   The  authors  felt  that  an  approximation 
of  this  type  which  goes  to  infinity  at  the  barrier  9  =  6-j 
would  fit  the  penalty  function  better  than  a  polynomial 
approximation  which  goes  to  infinity  only  as  9  ->°°.   9   is 
an  estimate  of  the  intersection  of  the  search  direction  s 
with  the  boundary  of  the  infeasible  region.   The  approximated 
penalty  function's  minimum  is  found  by  applying  Newton's 
method  to  find  a  step  length  which  corresponds  to  the  local 
minimum  along  the  direction  of  search. 

Since  the  explanation  of  this  method  by  the  authors 
was  rather  sketchy,  certain  parts  could  be  interpreted 
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differently  and  only  resolved  by  testing  the  method  on 
different  problems.   An  explanation  of  the  notation  for 
this  method  is  first  given. 


1   m+p     7  9 

F(8)  =  f(6)  +  ±   E    h/(9)  s  a  +  be  +  cBZ 


r  j=m+l  3 


m 
G(0)  =  -r   I   In  g. (0)  z 


j-1      3        (eu-0) 


sAVf (9)  +  £    Eh. (9)siVh. (9) 
F»(8)  =  r  j=m+l  3 J 


s  b  +  c9 


G*  (9)  =  - 


^_    m    sTvgi(9)    .       d 

s||  j-i   gj(9)     Oy-e)2 


The  single  variable  search  is  supplied  with  an  initial  point, 
functional  and  gradient  information  at  that  point  and  a 
search  direction.   The  steps  of  the  method  follow. 

STEP  1. 

A  feasible  step  length  in  the  direction  of  search  is 
needed  along  with  the  initial  point  in  order  to  estimate  the 
parameters  a,  b,  c,    d  and  0„.   A  unit  step  length  is  initially 
taken  in  trying  to  find  a  feasible  step  length.   If  it  is  not 
feasible,  the  step  length  is  progressively  halved  until  a 
feasible  step  length  is  found.   The  lower  and  upper  feasible 
step  lengths  are  written  as  9Q  and  9^,  respectively.   At  the 
start  of  the  search  9Q  =  0.   Function  and  gradient  evaluations 
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are  made  at  the  second  feasible  step  length  in  order  to 
compute  F^,    G,,  F  '  and  G,  '  (shortened  notation  for  F(61) , 
GO,)  ,  etc.)  . 

By  using  the  information  at  the  two  feasible  step 
lengths,  simultaneous  equations  are  solved  to  determine 
the  parameters .   The  parameter  values  are 


c  =  


2(9Q  -  e1) 


b  -  V  -  2c9o 


u  v  -v 


d  =  v(9u-eo)2  • 

The  parameter  6   is  taken  as  the  smallest  value  of  two 
computed  values  which  is  greater  than  6,. 

STEP  2. 

Once  the  parameters  are  determined,  the  equation  for 

the  minimum  of  T  which  is 

T1  =  0  =  b  +  2c9  +  j 

(eu-  er 

is  solved  by  Newton's  method  using  the  midpoint  between  the 
two  feasible  step  lengths  as  the  starting  point. 
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9l-90 
ei-0n       T'(   2   ) 


°1   °0 
T" (  x    u) 


2 

where  T"  =  2c  +  •=• 


3        3  * 

(eu"e) 


STEP  3 


The  values  of  VP(8)  are  computed  and  if  the  stopping 
criterion  is  satisfied,  the  single  variable  search  is 
terminated  with  9  as  the  optimal  step  length. 

STEP  4. 

If  the  stopping  criterion  is  not  satisfied,  a  new 
step  length  is  computed  using  Newton's  method  in  Step  2. 

The  authors  of  this  method  stated  that  different 
default  actions  would  be  necessary  in  cases  where  the  inter- 
polation failed  or  was  unsatisfactory.   For  example,  if  for 
the  feasible  step  lengths,  9fi  and  9,,  a  value  of  9   could 
not  be  determined,  it  would  be  necessary  to  use  another 
feasible  step  length  to  determine  where  the  intersection  of 
the  search  direction  and  the  boundary  of  the  infeasible 
region  was.   In  later  stages  of  the  minimization  the  authors 
said  that  this  method  sometimes  failed  and  a  quadratic 
interpolation  would  have  to  be  used. 

Implementing  this  method  in  the  SUMT  program  would 
also  require  more  storage  or  more  gradient  evaluations 
along  with  modifications  to  other  subroutines  and/or  additional 
subroutines. 
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No  comparisons  are  made  of  the  Lasdon,  Fox  and 
Ratner  or  Fletcher-McCann  methods  to  the  other  three 
methods  since  they  are  only  discussed  and  not  programmed 
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IV.   TEST  PROBLEMS 

A.  TEST  PROBLEM  ORIGIN 

The  problems  used  to  compare  the  three  different  single 
variable  search  methods  programmed  as  OPT  subroutines  were 
taken  from  problems  used  in  a  thesis  by  Lt.  J.  Waterman,  USN, 
which  compared  three  different  nonlinear  programming  codes 
[Ref .  7] .   The  SUMT  program  and  problem  data  decks  were  the 
same  as  those  used  by  Waterman  except  for  the  Golden  Section 
OPT  subroutine  being  replaced  by  the  cubic  and  quadratic 
interpolation  methods. 

The  structure  and  degree  of  difficulty  among  the  seven 
test  problems  are  quite  varied  and  represent  a  sample  of 
some  real  world  problems.   The  number  of  variables  and 
constraints  ranged  from  9  to  100  and  2  to  20  respectively. 
The  problems  contain  combinations  of  linear,  nonlinear, 
equality  and  inequality  constraints. 

The  following  subsection  contains  the  descriptions  of 
the  test  problems  as  presented  in  Waterman's  thesis. 
Constants  in  each  problem  are  represented  by  a,  b,  c,  d,  e,  L, 
m,  u  and  s  unless  otherwise  specified. 

B.  TEST  PROBLEM  DESCRIPTION 

Problem  1  was  an  example  of  determining  the  chemical 
composition  of  a  complex  mixture  under  conditions  of 
chemical  equilibrium.   It  contained  45  independent  variables 
and  16  linear  equality  constraints. 
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hk  x. 


minimize      f(x)    =     E       [      E      x.v(c.v   +    In  -r- — & — )     ] 

nk 
E       x 


'ikv    ik         ±n   ~h~ 
k=l        j  =  l      1*      1*  hk 


jk 
3=1      J 


hk 
7  K 


subject   to  h. (x)    =      E       (    E      E.  .,x.  ..  )    -  b.    =   0, 

k=l      j-1      ^k   ^k 

i   =   1, ... ,16 


xjk-°  j=l,...,hk        k  =   1, . . .,7 


Problem  2  was    formulated  by   the  Shell   Development 
Company.      It  consisted  of   15   variables   and   5   nonlinear  equality 
constraints . 


10  5        5  5  3 

maximize      f(x)=      E      b.x.-      E        E      yz-2      E      d.  z 

i=l      1   1        j=l   i=l  j=l      J 


where  y  =  cijx(1Q+i)   and   z  =  x(1Q+j) 


5  10 

subject  to  g.(x)  =  2   E   y  +  3  d . z  +  e .  -  E   a. .x.    0 

3         i=l         J  1  i=l  ±J    x 

x.  >_  0      i  -   1, . .  o  ,15 


Problem  3  was  to  maximize  the  area  of  a  hexagon  in  which 
the  maximum  diameter  was  unity.   The  problem  had  9  independent 
variables,  13  nonlinear  inequality  constraints  and  a  lower 
bound  of  0  for  xg . 

45 


Maximize:   f(x)  =  .5(XjX4  -  x2x3  +  x3xg  -  x5x9  + XgXg  -  XgX?) 


subject  to:      1  -  x2  -  x2  _>  0 

1  -  x42  >  0 
1  -  x5   -  xfi  >_   0 

1  -  {x1-x5)2   -  (x2-x6)2  >  0 

2  2 

1  -  (x1  -  x5)   -  (x2  -  xg)   _>  0 

2  2 

1  -  (x1-x?)   -  (x2-xQ)   >  0 

2  2 

1  -  (x3  -  x5)   -  (x4  -  x6)  >    0 

2  2 

1  -  (x3-x?)   -  (x4-Xg)   >  0 

2  2 

1  -  x?   -  Cx8  -  x9)   _>  0 

x,x4  -  x-x-  _>  0 

H 

-x5x9  >  0 

X5X8  "  X6X7  >  ° 

x9  >  0   . 


Problem  4  was  probably  the  most  difficult  of  the  test 
problems.   It  included  a  linear  objective  function,  24 
variables,  12  nonlinear  equality,  2  linear  equality  and 
6  non-linear  inequality  constraints.   The  variables  had 
to  also  be  zero  or  positive. 
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24 
minimize:    f(x)  =   E    a.x. 

i-1   x  X 


subject  to: 


h  (x)  =  (i+12) °iXi      =  0   i  =  1      12 

n^xj  24   x  12   x     o,  i   i,...,±^ 

b,.x1,,,   E     rJ-     40b.   E    r-i 

(i+12)    ,  _   b.        1-1   b. 

3=13    3  3=1        3 


24 

h,-(x)  -   E  x.  -  1  =  0 

13      i-1  x 

12  x.       24   x. 

hu(x)  =   E  -=i  +  f    E    ^  -  1.671  =  0 

±4      i=l  ai     i=13   i 


where  f  =  142.224 


- (x.  +  x t  .    , 0, ) 
h(i+14) <*>  "    *24   '     '   +  -i  >  °  .  i  -  1'2'3 

E   x. 
j-l  3 


m   -(X(i+3)  +  X(i+15))  +       o,i=4,5,6 
n(i+14)  ^x;  24  i  - 

E   x. 


x.  _>  0,  i  =  1,  ...  ,24  . 
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Problem  5  was  a  weapon  assignment  problem  with  100 
variables,  a  nonlinear  objective  function,  12  linear 
constraints  and  zero  lower  bounds  for  the  variables. 


20         5    x.  . 
minimize:     f(x)  =   S    u.  (   n    a.  3    -  1) 

j=l    3        i=l   ^ 


subject  to: 


I      x   -  b .  >  0     j  =  1,6,10,14,15,16,20 
i=l   XJ     J  _ 


20 

E   x.  .  +  c.    0     i  =  1, .  .  .,5 
j-1   1D 


xi-  _>  0      i  =  1,...,5,   j  =  1,...,20 


Problem  6  was  adapted  from  an  inventory  model  where 
the  x.,  i  =  1,...,50,  represent  the  reorder  quantity  for 
50  inventory  items  and  x.,  i  =  51,..., 100,  represent  the 
reorder  points  for  the  same  50  items.   It  contained  100 
variables,  1  linear  and  1  nonlinear  inequality  constraint, 
and  50  lower  bounds  on  the  variables. 
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50    B.  (x.  +  50) 
minimize:     f(x)  =   £    — 1 


i=i         xi 


where 


1     n  o    <*.     s .d.    d. 

B.(xi  +  50)  =  i  (s.2*^2)  #(-1)  -  -Jj-i  ♦(ji)  , 


1 


2 

di  =  x(i+50)  "  mi  '  ♦■(*)  ■  7=-  e"x  /2 


'2tt 


X 

and   $(r)  ■  /   (J)(x)  dx  . 

t 


subject  to: 


50      x. 

200,000  -   Z    c.(^+  x(.+50)  -  m.)  >  0 


50  L. 

300  -   Z  —  >  0 

.  ,  x.  — 
1=1    l 


xi  >  0    ,      i  =  1, . . . ,50. 


Problem  7  was  an  entropy  model.   There  were  4  6  population 
centers  connected  by  a  transportation  network.   Using  a 
congestion  cost  function  the  model  yields  an  equilibrium 
solution  that  identifies  nodal  populations  as  entropic 
functions  of  the  total  cost  of  the  journey  to  work.   The 
problem  contained  4  6  variables,  all  of  which  have  a  zero 
lower  bound,  1  nonlinear  inequality  and  1  linear  equality 
constraint. 
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minimize:    *<*>-*    SOT  <*»  JOT  > 


46 
subject  to:      500  -   E   x.  =  0 

i=l   1 


46  h 

10000  -   E   c.y.  +  ad.y.   >  0 

i-1   X  x     x  l  " 


where     y.  =  x.  +    I  x.    ,    A(i)  consists  of  all  the 

jeA(i)  3 

arcs  that  converge  directly  and  indirectly  upon  node  i, 

and 

x^>^0      i  =  l,...,46. 


The  preceding  problem  descriptions  are  only  meant  to 
give  a  feeling  of  the  kind  of  test  problems  used  and  their 
structure.   Detailed  information  of  how  each  problem  was 
set  up,  the  constant  values,  initial  starting  points  and 
where  the  problems  specifically  originated  is  contained 
in  Ref.  5. 
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V.   TEST  RESULTS  AND  COMPARISONS 

A.   PROGRAMMING  AND  TESTING  PROCEDURE 

The  cubic  interpolation  method  was  first  programmed 
as  a  program  that  did  a  single  variable  search  when  given 
a  penalty  function,  its  gradient,  an  initial  starting  point 
and  a  search  direction.   After  it  had  been  debugged  it  was 
converted  into  an  OPT  subroutine.   After  the  subroutine 
was  debugged  using  a  couple  of  the  less  difficult  test 
problems,  OPT  was  converted  into  the  quadratic  interpolation 
method  by  modifying  the  bracketing  procedure  and  changing 
the  interpolation  from  a  cubic  function  to  a  quadratic 
function. 

When  the  quadratic  interpolation  method  had  been  debugged, 
each  problem  was  run  through  the  SUMT  program  three  times, 
using  the  three  different  search  methods.   As  discussed 
earlier,  the  stopping  criteria  was  made  essentially  the  same 
in  each  of  the  OPT  subroutines  so  that  the  only  thing  that 
varied  in  solving  each  test  problem  was  the  single  variable 
search  method. 

In  the  process  of  running  the  other  test  problems,  minor 
debugging  changes  had  to  be  made  in  the  interpolation  methods. 
After  all  the  changes  had  been  incorporated,  a  final  run 
of  each  problem  with  each  method  was  made  which  produced  the 
results  used  in  the  comparisons.   Each  method  found  the  same 
solution  to  each  problem  to  within  six  significant  figures. 


51 


B.   COMPARISON  CRITERIA 

Perhaps  the  best  way  to  compare  the  different  search 
methods  is  to  compare  the  computation  time  required  to  find 
the  solution  to  the  problem.   However,  because  of  computer 

interactions,  computation  time  may  vary  as  much  as  twenty 

> 
five  per  cent  when  the  same  problem  is  run  at  two  different 

times .   Computation  times  were  computed  for  each  problem  to 

see  if  a  trend  could  be  seen  between  the  methods. 

Counters  were  inserted  in  the  SUMT  program  so  that  the 
number  of  single  variable  searches  performed  and  the  number 
of  functional  and  gradient  evaluations  needed  to  find  the 
solution  could  be  counted. 

Results  of  the  test  problem  runs  are  shown  in  Table  I. 
In  running  problem  7  with  the  interpolation  methods  it  was 
necessary  to  change  the  factor  by  which  the  increment  was 
multiplied  for  step  length  increase  or  reduction  to  1.618 
vice  2  and  0.618  vice  0.5  respectively.   This  change  was 
needed  because  a  feasible  starting  point  was  found  in  the 
bracketing  procedures  which  the  SUMT  program  could  not  solve 
By  changing  the  factors  to  the  same  as  those  in  the  Golden 
Section  method,  the  same  feasible  starting  point  was  used 
to  solve  the  problem. 

C.   COMPARISONS 

Theoretically,  each  single  variable  search  should  reach 
the  same  solution  for  all  three  methods.  However,  because 
of  the  tolerances  allowed  in  the  stopping  critera,  the 
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TABLE  I.   TEST  PROBLEM  RESULTS 


Problem 


TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 

TIME 
SVS 
F 
G 


TIME  =  Computation  time 
SVS  =  #  single  variable  searches 
F   =  #  function  evaluations 
G   =  #  gradient  evaluations 


Golden  Section 

Cubic 

Quadratic 

105.8 

105.7 

94.8 

86 

89 

75 

15487 

6171 

10149 

1462 

7497 

1275 

11.4 

9.6 

8.3 

63 

53 

53 

4354 

1574 

2606 

378 

1704 

318 

4.8 

5.6 

4.0 

39 

38 

40 

8606 

3302 

5862 

613 

3358 

628 

395.3 

404.1 

396.1 

163 

151 

174 

44255 

13922 

30048 

92.27 

97319 

98694 

103.7 

97.9 

97.0 

20 

19 

19 

2683 

1250 

1991 

296 

1415 

271 

165.4 

145.9 

155.5 

36 

34 

35 

2320 

836 

978 

104 

817 

107 

78.8 

80.3 

67.8 

20 

20 

17 

865 

307 

397 

1993 

2257 

1708 
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solutions  were  slightly  different  at  the  conclusion  of  each 
single  variable  search.   This  meant  that  at  the  start  of 
the  next  single  variable  search  the  search  direction  would 
be  slightly  different  also.   As  several  single  variable 
searches  were  performed  in  each  problem,  the  differences 
accumulated  and  this  explains  why  for  some  of  the  problems 
it  takes  a  different  number  of  searches  for  each  method  to 
reach  the  same  solution.   This  also  explains  the  difference 
in  the  number  of  function  and  gradient  evaluations  required 
for  each  problem  for  the  Golden  Section  and  quadratic 
interpolation  methods. 

Table  II  shows  the  normalized  results  for  the  time  per 
single  variable  search,  function  evaluations  per  single 
variable  search  and  gradient  evaluations  per  single  variable 

As  was  expected,  the  number  of  gradient  evaluations  per 
single  variable  search  is  essentially  the  same  for  Golden 
Section  and  quadratic  interpolation -methods.   The  only  thing 
that  can  be  compared  between  the  three  methods  is  the 
computation  time  per  single  variable  search.   Function 
evaluations  per  single  variable  search  can  only  be  used  as 
a  measure  of  effectiveness  for  Golden  Section  and  quadratic 
interpolation  since  the  cubic  interpolation  also  requires 
gradient  evaluations  and  fewer  function  evaluations. 

In  looking  at  the  times  per  single  variable  search  the 
quadratic  interpolation  was  faster  than  the  Golden  Section 
and  cubic  interpolation  methods  5  out  of  7  and  4  out  of  7 
test  problems  respectively.   Cubic  interpolation  was  faster 
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TABLE  II.   NORMALIZED  TEST  PROBLEM  RESULTS 


Problem 


Golden  Sectioi 

l    Cubic 
1.19 

Quadratic 

1.23 

1.26 

180.1 

69.3 

135.3 

17 

84.2 

17 

.181 

.181 

.157 

69.1 

29.7 

49.2 

6 

32.1 

6 

.123 

.147 

.1 

220.7 

86.9 

146.6 

15.7 

88.4 

15.7 

2.42 

2.67 

2.28 

271.5 

92.2 

172.7 

565.2 

644.5 

567.2 

5.18 

5.15 

5.10 

64.4 

24.6 

27.9 

2.8 

24.0 

3.0 

4.59 

4.29 

4.44 

64.4 

24.6 

27.9 

2.8 

24.0 

3.0 

3.94 

4.02 

3.98 

43.3 

15.4 

23.4 

99.7 

112.9 

100.4 

TIME 

F 

G 

TIME 

F 

G 

TIME 

F 

G 

TIME 

F 

G 

TIME 

F 

G 

TIME 

F 

G 

TIME 

F 

G 


TIME  =  time/single  variable  search 

F    =  #  function  evaluations/single  variable  search 

G    =  #  gradient  evaluations/single  variable  search 
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than  Golden  Section  in  3  of  the  7  problems  with  the  same 
time  for  problem  2 . 

In  comparing  the  number  of  function  evaluations  per 
single  variable  search,  the  quadratic  interpolation  method 
required  fewer  than  the  Golden  Section  on  all  seven  problems. 
In  most  of  the  test  problems,  it  required  between  thirty  to 
forty  per  cent  fewer  function  evaluations.   The  cubic  inter- 
polation method  required  fewer  function  evaluations  than  the 
other  two  methods  because  it  used  more  information  about  the 
penalty  function  at  each  feasible  step  length  to  find  the 
optimal  step  length.   However,  it  required  more  gradient 
evaluations  which  increased  the  computation  time  making  it 
about  the  same  as  the  Golden  Section  and  slightly  slower  than 
the  quadratic  interpolation. 

With  a  larger  sample  of  test  problems  more  statistically 
sound  comparisons  could  be  made  of  the  three  different  methods 
Also  samples  of  test  problems  that  are  similar  in  structure 
could  be  tested  to  show  if  one  method  worked  better  than  the 
others  on  those  type  of  problems. 

It  should  be  emphasized  that  the  results  in  this  thesis 
are  obtained  from  single  variable  searches  on  a  very  special 
type  of  function  —  the  unconstrained  penalty  function  in  the 
mixed  interior-exterior  penalty  function  algorithm.   No 
attempt  should  be  made  to  generalize  these  results  to  other 
nonlinear  programming  methods  which  also  use  single  variable 
searches,  but  on  a  significantly  different  class  of  functions. 
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VI.   CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  STUDY 

A.   CONCLUSIONS 

In  looking  at  the  computation  time  required  for  each 
single  variable  search,  the  quadratic  interpolation  method 
was  faster  than  the  other  two  methods  on  the  test  problems 
that  were  run.   The  quadratic  interpolation  method  uses 
information  about  the  penalty  function  at  each  feasible  step 
length  to  reduce  the  interval  of  uncertainty  whereas  the 
Golden  Section  reduces  the  interval  of  uncertainty  by  a 
constant  factor.   Because  of  this  difference,  quadratic 
interpolation  proved  to  be  slightly  more  efficient  than 
Golden  Section.   The  only  drawback  in  the  quadratic  inter- 
polation method  is  the  fact  that  bracketing  the  minimum  can 
be  quite  hard  to  do.   The  cubic  interpolation  method  reduced 
the  interval  of  uncertainty  in  a  more  efficient  manner  than 
quadratic  interpolation  or  Golden  Section.   However,  having 
to  make  gradient  evaluations  to  determine  the  directional 
derivative  values  proved  to  be  very  costly.   The  time  per 
single  variable  search  was  about  the  same  as  the  Golden 
Section  search  method. 

If  the  user  of  the  single  variable  search  method  has 
prior  knowledge  about  the  complexity  of  the  function  or 
gradient  evaluations  he  may  prefer  cubic  interpolation  over 
Golden  Section  or  quadratic  interpolation,  or  vice  versa, 
since  the  cubic  interpolation  requires  more  gradient  evalua- 
tions and  fewer  function  evaluations. 
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B.   SUGGESTIONS  FOR  FURTHER  STUDY 

A  combination  of  different  methods  made  into  a  single 
variable  search  could  be  done  with  the  three  different 
methods  that  were  programmed.   This  combined  methods  approach 
could  use  Golden  Section  until  the  minimum  was  bracketed, 
then  use  quadratic  interpolation  until  the  minimum 
was  bracketed  much  closer  and  then  finally  use  the  cubic 
interpolation  to  home  in  on  the  local  minimum  of  the  search. 
The  Lasdon,  Fox  and  Ratner  method  uses  three  different  stages 
that  each  find  a  feasible  step  length,  with  the  final  stage 
finding  an  optimal  step  length  for  the  search.   Another 
approximation  of  the  penalty  function  like  the  Fletcher- 
McCann  method  could  also  be  devised. 

Implementation  of  the  Lasdon,  Fox  and  Ratner  or  Fletcher- 
McCann  method  into  the  SUMT  program  with  runs  of  the  test 
problems  would  enable  a  comparison  to  be  made  not  only  to 
the  three  methods  tested  in  this  paper  but  also  between  the 
Lasdon,  Fox  and  Ratner  method  and  Fletcher-McCann  method. 

Another  possibility  for  study  would  be  to  fine  tune  the 
three  different  methods  by  adjusting  the  tolerances  and 
termination  conditions  to  see  that  if  by  allowing  a  less 
accurate  determination  of  the  minimum  of  the  single  variable 
search,  the  solution  to  the  unconstrained  problem  could  be 
found  any  faster. 

If  a  larger  sample  of  test  problems  including  more  and 
varied  types  of  nonlinear  programming  problems  could  be 
tested  using  the  three  different  methods,  a  better  comparison 
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could  be  made.   The  results  may  either  show  one  method  to 
be  superior  over  the  other  methods  in  all  of  the  problems 
or  show  each  method  suited  best  to  a  specific  type  of 
problem  enabling  the  user  to  choose  the  method  which  best 
applies  to  the  situation  at  hand. 
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APPENDIX 
GLOSSARY  FOR  GOLDEN  SECTION  OPT  SUBROUTINE  [Ref.  2] 


DELX 

DELXO 

DOTT 

I 
ISW 

J 

KSW 

M 

MN 

N 

NSATIS 

NTCTR 

N401 

N404 

N405 

PREV3 


The  vector  indicating  the  direction  of  move 
in  one  dimensional  optimization.   S 

The  gradient  vector  of  the  penalty  function 
VP. 

The  inner  product  of  the  move  vector  and  the 
negative  gradient  vector  of  the  penalty 
function.   -sTVp 

Used  as  an  index  in  DO  loops. 

A  switch  showing  whether  the  motion  on  a  given 
vector  failed  to  decrease  the  penalty  function 
and  the  negative  was  tried. 

Used  as  an  index  in  DO  loops . 

Switch  showing  whether  less  than  25  feasible 
points  were  found  on  the  direction  vector. 

The  number  of  inequality  constraints 

Number  of  moves  in  search  for  a  solution  of 
a  subproblem 

The  number  of  variables 

Indicates  whether  constraints  are  satisfied. 

The  number  of  the  point  on  which  the  program 
is  working. 

The  number  of  points  generated  in  attempting 
to  find  a  point  along  the  search  direction. 

The  number  of  infeasible  points  found  on  the 
direction  vector. 

Switch  showing  that  a  feasible  point  on  the 
direction  vector  could  not  be  found  and  the 
negative  was  tried. 

Penalty  function  value  at  lower  step  length 
P(X3) . 
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PX1  Penalty  function  value  for  left  interior  steplength 
P(X2). 

PO  Current  value  of  the  penalty  function   P(X) 

PI  Penalty  function  value  for  upper  steplength 
P(X1) 

P31  Penalty  function  value  at  initial  starting  point 

RJ  Vector  of  current  values  of  the  constraints 

RJ1  Vector  of  previous  value  of  the  constraints 

X  Vector  of  current  value  of  the  variables   xq+^s 

XX  Temporary  storage  used  for  switching  vector  values 

XI  X  vector  of  upper  step  length  X_+0us 

X2  X  vector  of  left  interior  step  length  XQ+6as 

X3  X  vector  of  lower  step  length  XQ+0Ls 
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SLERCUTINE  OPT 

INFLICIT  REAL*8(A-H,0-Z) 

C  *ARCH  1971 

C 

C  0P7  LCCKS  FCR  A  MINIMUM  ALCNG  THE  SEARCH  VECTOR  USING  THE 

C  GCLDEN  SECTION  SEARCH  METHOD 

CCMMCN/SHARE/  X(IOO),  CELUOO),  A  ( 100, 10G  )  i N , H ,  MN,NF1 

f£EK££!,^S!sY£/  F,G,PO.RSIGMA,  RJ(2G0),  RHC 

CCNMCN/CRST/  DtLX(lOO),  DELXO(IOO),  RH0IN,R£7I0,  EPSI , 
1THETAO, 

2  RSIG1,  Git     Xl(lOO),  X2(100),  X2(100),  XR2Q0O), 
3XP1 (100 ) tPRlt 

£  Fg2, PI j  Fit  RJl(200)t  DOTT,  PGRAC(IOO),  CIAG(IOO), 
5  FPEV3,ADELX,  NTCTR,  NUMINI,  NPHASE,  NSA7IS 

AES(CUMMY)=DABS(OUMNY) 

KSM1 

N AC5=1 

F31=P0 

I  SV%  =  1 

CCTT=0. 

DC  10  J  =  1,N 
10     CCTT=C.7T+DELX( J)*DELXO(J) 

GC  TO  40 
20     CC  20  I  =  1,N 
20     CELX(I) =-OELX(i) 
40     CONTINUE 

,N^C4=0 

*<N  =  KN  +  1 
C  MN  IS  NOW  NUMB*  OF  P0IN7S  AFTER  MIN  ACFIEVEC 

NTCTR=NTCTR+1 

CC  50  1  =  1, N 
50     X2(I)=X(ll 

FX1=P0 

N401=0 
60  MC1  =  N4C1  +  1 

CC    70    I=1,N 
70  X<I)=X2<I )+DELX(I ) 

CALL  EVALU 
C  1  VEAIVS  SA7IS.GF  CONSTRAINT  NT.PREV.  2MEANS  NCCFANGE  2 
C  MEANS  VIOLATION 
C  IF  POINT  IS  NOT  FEASIBLE  GIVE  IT  AN  ARBITRARILY  HIGH  VALU 

GC  TO  (540,90,80),  NSATIS 
60     P>2=10.E25 

FO10.E25 

GC  TO  100 
90     CONTINUE 

FX2=P0 

IF  (PX1-PX2)  100,100,150 
100    IF  (N401-2)  130,110,110 
110    CC  120  1=1, N 

12C      xim=xm 

F1=PX2 

GC  TO  420 
C  ONLY  CNE   PCINT   SO  FAR  COMPUTEC 
120    CC  140  1  =  1, N 
140    X2(I)=X2(I) 

FREV2=P>1 

GC  TO  160 
150    CC  160  1  =  1, N 

X2(I)=X2(  I) 

K*(I  I"KII  ) 

160    CELX(I)=1.61803399*CELX(I) 

PFEV2=PX1 

F)1=PX2 

GC  TO  60 
C  GCLOEN  SECTION  SEARCH  METHOD. 
C  E  VECTCR   GCES  TO  XKI  ) 
170    PC=1.E2  6 

MC4  =  N404+1 
180    CC  190  1=1, N 


190    X1(I)=X(I) 

F1=F0 

CC  200  1=1, N 

X(I)=.28196601*(X1(I)-X2( I) )+X3(I) 
200    X2(I>=X(I) 

CALL  EVALU 

GC  TC  (540,270,210),  NSATIS 

210  IF  (N404.LT. 30)  GC  TC  170 

211  CCNTINUE 

C  THERE  IS  NG  REFERENCE  TC  211,  TH5  AeCVE  STATEMENT  IS  A 

C  DUMMY  STATEMENT 

C—  IT  IS  POSSIBLE  NG  FEASIBLE  POINT  EXIST,  IF  NCT  TRY  MGVING 

C  ON  CELXO 

C —   IF  IT  IS  NCT  POSSIBLE  TC  MOVE  GN  CELXC  THEN  kc  MUST  EE 

C   AT  A  SCLUTION  CF  NIP  PROBLEM. 

IF  (N404.GT.100)  GO  TO  240 
220    CC  230  1  =  1, N 

IF  (AeS(ABS(X3(I  )/XK  I)  )-l.  ).GT.l.E-7)  GC  TC  170 
220    CCNTINUE 
240    GC  TC  (25C260),  N405 
250    N4C5=2 
C — TRY  TC  MCVE  CN  GRADIENT. 

NTCTR=NTCTR-1 

m=NN-l 

GC  TC  20 
260    WRITE  (6,580) 

CALL  TIMEC 

CALL  OUTPUT  (1) 

CALL  REJECT 

STCP  22042 
C 
270    CCNTINUE 

N*C4=0 

FX1=P0 

CC  280  1  =  1, N 
280    MI  )=G.2319  6601*(X1(I  )-X2(I))+X2(I) 

CALL  EVALU 

GC  TO  (540,290,220),  NSATIS 
290    FX2=P0 

N401=l 
200    N^C1=N401+1 

IF  (N401-25)  340,210,210 
210    KSW  =  2 

IF  (N4C1-4Q)  320,460,460 

IF  UBS(X2(I)/X(I)-1.0).GE. 1.5-7)  GO  TO  240 
220    CCNTINUE 

GC  TC  4  60 
240    IF  (ABS(PX1/PX2-1. ULE.l.E-7)  GC  TC  460 

IF  (FX1-PX2)  350,46C,400 
C  FRCM  IEFTGRIGHT   X3 (  I  )  < PREV2 ) X2 ( I  )  ( PX  1  )X ( I ) PX2  X  1  ( I  )  P  1 
250    CC  260  1  =  1, N 
260    X1(I)=X(I) 
C  TFRCW  A*AY  RIGHT  PART 

F1=PX2 

CC  270  1  =  1, N 
C  PCINTXF1  E5CCMES  XP2  „  , 

270    X(I)=.3  81966G1*(X1(I)-X3(I))+X3(!) 
C  TEMPORARILY  IN   X  STORAGE 

GCLT0E^(540,280,170),  NSATIS 
280    CCNTINUE 

c  sutchSectcrs  to  proper   POSITION 

PX1=P0 

CC  290  I=1,N 

XX  =  X2(I  ) 

X2(I)=X(I) 
290    X(I)=XX 

GC  TC  200 
C  Li=FT  SICE   TOSSED   AWAY 
c — "CHANGFS  FCR  NONUNIMCCAL  FN 
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c— 

400 
410 

420 


420 
44C 


450 


C  T 
460 


HE 


470 


480 
490 
500 

C  IF 

510 
5  20 
520 


C  A 
540 

550 


C  — 

c  - 

560 

570 

C 
5  80 


GC  7C  ThRCW  AWAY  RIGHT  IN  THIS  CASEINIT  VAL  LT  FI3  PT 
IF  (PREV2-PX2)  35C25C410 
CC  420  I-1,N 
X2(Ii=X2(I  ) 
X2(I)»X(I) 
FPEV3=PX1 
F)1-PX2 
CC  440  1  =  1, N 

X(I  )=0.281S6601*(X1(I *-X2(I))+X2(I) 
CALL  EVALU 

GC  TC  (540,450,170),  NSATIS 
CCNTINUE 
FX2=P0 
GC  TC  300 

UTERICR  PGINTS  NOW  GIVE  EQLAL  VALUE  FOR  P.  CCMPUTE  MI 
CC  470  1  =  1, N 
CELXOd  )  =  X(I  ) 
X(I)  =  (CELXC(I  )+X2(I) )*0.5 
CCNTINUE 
CALL  EVALU 

GC  TO  (460,490),  KSW 

IF  (ABS(P0/PXl-U).GT.lBE-7)  GO  TC  520 
GC  TO  (500,510),  ISI« 
IF  (P0.LT.P21)  GO  TO  510 
ISk=2 
P-FLNCTICN  CIDN,T   GC  DOWN  TRY  NEG  VECT. 
GC  TO  20 
RE7LFN 

CC  530  1  =  1, N 
X(I)=OELXO(I) 
GC  TO  250 

VE    NOW  IN  FEASIBILITY  PHASE 
CC  550  1=1, M 
IF  (RJ(I))  560,5oO,550 
CCN7INLE 
NSATIS=4 
PETLPN 
PPCELEM  HAS  BECOME  FEASIBLE 

F  -  FUNCTION  CHANGES  IF  A  CONSTRAINT  BECC.VFS  FEASIELE 
f*h'Q 

CC  570  1=1, M 
R„1(I)=PJ(I) 
RETLRN 

FCFNAT  (  80H   OPT  CAN-T  FINC  A  FEASIBLE  POINT, THAT 
1GIVES  A  LOWER  VALUE  OF  THE  F-FUNCTION.  ) 

ENC 


PE 
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GLOSSARY  FOR  CUBIC  INTERPOLATION  OPT  SUBROUTINE 

ADELX    Magnitude  of  the  search  direction  vector 

AL       Lower  step  length  Q_ 

BL       Upper  step  length   0 

COTES    Cosine  of  the  angle  between  the  gradient  and  search 
direction  vectors 

D        Factor  by  which  increment  is  multiplied 

DELX     Search  direction  vector   s 

DELXO    Gradient  vector  of  the  penalty  function   VP(X) 

DLX1     Gradient  vector  of  penalty  function  for  lower  step 
length   VP(X3) 

DLX2     Gradient  vector  of  penalty  function  for  upper  step 
length   VP(X2) 

DOTS     Directional  derivative  of  interpolated  step  length 
times  ADELX  P1 (X) 

DOTT     Directional  derivative  of  initial  starting  point 
times  ADELX   P1 (XI) 

DOTTA    Directional  derivative  of  lower  step  length  times 
ADELX  P' (X3) 

DOTTB    Directional  derivative  of  upper  step  length  times 
ADELX   P' (X2) 

EPSO  Stopping  criteria  tolerance  for  cosine  test 

I  Index  used  in  DO  loops 

J  Index  used  in  DO  loops 

K  Index  used  in  DO  loops 

KTER  Number  of  points  evaluated  in  bracketing  the  minimum. 

KTR  Number  of  cubic  interpolations  performed 

M  Number  of  inequality  constraints 

MN       Number  of  moves  in  search  for  solution  of  a 
subproblem 
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N  Number  of  variables 

NRED  Number  of  consecutive  step  length  reductions 

NSATIS  Indicates  whether  constraints  are  satisfied 

NTCTR  Number  of  the  point  on  which  the  program  is  working 

N405  Switch  showing  that  a  feasible  point  on  the 
direction  vector  could  not  be  found  and  the 
negative  was  tried. 

PO  Penalty  function  value  at  current  xvector  P(X) 

PX1  Penalty  function  value  at  lower  step  length  P(X3) 

PX2  Penalty  function  value  at  upper  step  length  P(X2) 

Q  Coefficient  used  in  computing  interpolated  step  length 

RJ  Vector  of  current  values  of  the  constraints 

RJl  Vector  of  previous  values  of  the  constraints 

SMAG  Magnitude  of  the  gradient  of  the  penalty  function 

TAL  Current  step  length  0 

TINC     Increment  by  which  step's  length  increased  or 
decreased 

X  Current  X  vector 

XI  Initial  starting  point  X  vector 

X2       X  vector  of  upper  step  length   XQ+0uS 

X3       X  vector  of  lower  step  length   Xo+0LS 

Z        Coefficient  used  in  computing  interpolated  step 
length 
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„    m         SIEFCUTINE    OPT 

r    J^TuF£TA?y?£Cy7£NETUSES    CUBIC    INTERPOLATION    TC    FIND    THE 
C    flMHf    ALONG    THE    GIVEN    SEARCH    VECTOR 

IMPLICIT  RCALS'C8(4—  H  1  —  Z) 

CCMMCN/SHARE/  X(IOO)1,  CEL(IOO),  A  ( 100,  100  )  t  N,M,  NN,NP1 

^CMNON  /VALUE/  F , G ,PQ .RS IGMA,  PJ(2C0),  PHC 
CCMNQN/CRbT/  DELXQOO),  DELXO(IQO),  RHOIN, RATIO,  EFSI, 
1< FETAO, 

3XP1  (100  )GPR1    X1(10C,T  x2dO0),  X2Q00),  XR2(100), 

i    SSi.'^lI^lt'RJKZOOt  DOTT,  PGRAC(IOO),  CIAGUOO), 
5.?Piyi»ADfcLX»  NTCTR,  NUMINI»  NPHASE,  NSATIS 

CINENSICN  CLX1(100),CLX2(10C> 

AES(CUMMY)=CABS(DUNMY) 

S C PT ( DUMMY )=DSQRT( DUMMY) 

EFSC=l„E-7 

CCTS=O.C 

N4C5=1 

GC    TO    5 
2    CC    4    1=1, N 
A    CELX(I)=-DELX(I  ) 
5    CONTINUE 
C    INITIALIZATION 

NMNN  +  1 

NTCTR=NTCTR+1 

CCTTA=0.0 

TAL=0.0 

AL=C.O 

C  =  2.0 

KTER=0 

NFEC=0 

TINC=1.0 

SMC-=O.C 

CC  10  K«1,N 

Sf AG=SMAG+CELX(K)**2 

CCTTA=CCTTA-DELX(K)*DELXO(K) 

XKK)-X(K) 

CLX1(K)=DELX0(K) 
1C  X2(K)=X(K) 

CCTT=-CCTTA 

FX1=P0 
C  INITIAL  ERACKETING  BEGINS 
25  TAL=TAL*TINC 
2C  CC  25  K  =  1,IS 
25  X(K)=X1(K)+DELX(K)*TAL 

KTER=KTEP+1 

CALL  6VALU 

GC  TO  ( 170,50,40)  ,  NSATIS 
C  RECLCE  STEP  SIZE 
4C  C=C5 

NFEC=NRED+1 

TINC=(TAL-AL)*D 

T/L=TAl-TINC 

GC  TO  30 
C  SECCNC  FEASI8LE  POINT  FCUND 
5C  Bl=TAL 

FX2=P0 

CALL  GRAD(2) 

CCTTB=O.C 

CC  60  K=ltN 

X2(K)=X(K) 

CLX2(K)=DSLX0(K) 
tC    CCTTE  =  CCTT6-DELX(K)*DELX0(K  ) 
C  CHECK  STOPPING  CRITERIA  OR  MAYBE  REVERSE  SEARCH  DIRECTION 

I  MAOELX.EC.0.0)  GC  TO  91 

IF(SMAG.EC.O.O)  GC  TO  91 

CCTES=  ABS(DOTTB)/(SQRT<SMAG)*ACELX) 

IF  (COTES. LT.EPSQ)  GO  TO  91 

IF(KTER.GT.IOO)  GO  TC  280 

IF(NREC.GT.20)  GO  TO  280 

CC  75  K  =  1,N 
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IF(ABS(ABS(Xl(K)/X2(K))-l.).GT.l.E-7)  GO  TC  76 

GC  TO  280 
_   76  IF(AES<A6S<PX2/PXl)-l.).LE.l.E-7)  GO  TO  91 
C  IF  CEPIVATIVE  IS  POSITIVE  THEN  THE  MINIMUM  13  8FACK3TEC 
C  ANC  CLEIC  INTERPOLATION  CAN  BE  PERFORMED,  IF  NOT  NAKS  THIS 
C  PCINT  THE  LCWER  STEP  LENGTH  AND  STEP  CLT  FARTHER 

IF(DOTTE.GT.O.O)  GQ  TC  100 

NFEC=0 

TINC  =  TIfsC*D 

CC  60  K=1,N 

DLXl(K)rDLX2(K) 
EC  X3(K)=X2(K) 

F>1=PX2 

Al^EL 

CCTTA=CCTTB 

GC  TO  2  5 
C  OF  THE  TWO  FEASIBLE  POINTS  RETURN  WITH  THE  SMALLEST 
SI  IF(FX2.LE,FX1)   GC  TO  220 

PC=FX1 

CC  92  K=1.N 

X(K)=X3 (K) 

CELX0(K)=DLX1(K) 
92  CONTINUE 

GC    TO    220 
C    CLEIC    INTERPOLATION    FOLLOWS    ************************* ***** 
ICC    KTP=0 

ics  ccminue 

Z=3.*<(  FX1-PX2)/(EL-AL) )+DOTTA+CCTTB 

C=SCRT(Z**2-DOTTA*CCTT8) 

TAL=BL-(BL-AL)*(DCTTB+Q-Z)/(DGTTB-CCTTA+2.*C) 

CC  130  K  =  1,N 
13C  X  (K)  =  X1  (K)+DELX(K)*T/SL 

CALL  EVALU 

GC  TO  (170,135,40) ,  NSATIS 
135  CALL  GRAD12) 

CCTS=O.C 

CC  140  K=1,N 

CCT5=DC7S-CELX(K)*0SLX0(K) 
14G  CONTINUE 
C  CHECK  STOPPING  CRITERION  CR  IF  SEARCH  CIRECTICN  SHCULC  EE 
C  REVEPSEC 

I F(ADELX.EC.O.O)  GC  TO  471 

IF(SMAG.EC.O.O)  GC  TO  471 

CCTES=AES(CCTS)/(SCR7(SMAG)*ADELX) 

IF  (CGTES.LT.EPSO)  GO  TC  471 

IF(ABS(ABS(PX2/O0  )-l. ) ,L E. l.E-7  )  GO  TO  471 

K7R*KTR+1 

IF (KTR.GT.10)  GO  TO  471 

CC  143  K=1,N 

IF(ABS(ABS(X2(K)/X(K)  )-l. ) .G7. 1. E-7 )  GO  7G  163 
143  CONTINUE 

GC  TC  220 
163  IF(COTS.GT.O.O)  GO  TO  150 
C  TFPCfe  AWAY  LEFT  SIDE 

AL=TAL 

FX1=PC 

CCTTA=CCTS 

CC  145  K=1,N 

CLX1(K)=DELX0(K) 
145  X3(K)=X(K) 

GC    TO  105 
C  THFCW  AWAY  PIGHT  SIDE 
15C  EL=TAL 

FX2=P0 

CCTT8=CCTS 

CC  160  K=1,N_ 

CLX2(K)=DELX0(K) 
16C  X2(K)=X(K) 

GC  TO  105 
170  CC  130  K  =  1,M 

IF(RJ(KJ)  190,190,160 
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160  CCNTINLE 

N5ATIS=4 

RETLRN 
ISC  PN  =  C 

CC  2C0  K=1,M 
2CC  RJ1(K)=PJ(K) 

PETLRN 
RETLRN  luITH  SMALLEST  VALUE 

471  IF(FO.LE.PXl)  GO  TO  474 
IF(FX1.LE.PX2)  GO  TO  475 

472  FC=FX2 

CC  473  K  =  1,N 
X(K)=X2<K) 
CELXO(K )=DLX2(K) 
472  CCNTINUE 
GC  TO  2  20 

474  IF(F0.LE.PX2)  GO  TC  220 
GC  TC  472 

475  PC=FX1 

CC  476  K  =  1,N 

X(K)=X2(K) 

CELX0(K)=DLX1(K) 

476  CCNTINUE 
22C  CCNTINUE 
225  CCNTINUE 

RETLRN 
REVERSE  SEARCH  DIRECTIONS 
26C  GC    TC  ( 290,200)»N405 
29C  N4C5=2 

NTCTR=NTCTR-1 

fMNN-1 

GC  TO  2 
2CC  kFITE<6»580) 

C4LL  TIMEC 

C*LL  OUTPL'T(l) 

CALL  REJECT 
5£0  FCFMAT(  80H   OPT  CAN-T 
1G1VES  A  LOWER  VALUE  OF 

STCF    22042 

EISC 


FINO  A  FEASIBLE 
THE  F-FLNC7IGN. 


PGINTfTHAT 
) 
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GLOSSARY  FOR  QUADRATIC  INTERPOLATION  OPT  SUBROUTINE 

AL       Interpolated  step  length 

D        Factor  by  which  increment  is  multiplied 

DELX     Vector  indicating  the  direction  of  move  in  one 
dimensional  optimization   S 

DELXO    Gradient  vector  of  the  penalty  function   VP 

DOTT     Directional  derivative  of  initial  starting  point 
times  search  direction  magnitude   P'(X1) 

I        Index  used  in  DO  loops 

J        Index  used  in  DO  loops 

K        Index  used  in  DO  loops 

KTER     Number  of  points  evaluated  in  bracketing  the 
minimum 

KTR      Number  of  quadratic  interpolations  performed 

M        Number  of  inequality  constraints 

MN       Number  of  moves  in  search  for  solution  of  a 
subproblem 

N        Number  of  variables 

NFIRS    Indicates  if  more  than  two  feasible  step  lengths 
found 

NRED     Number  of  consecutive  step  length  reductions 

NSATIS    Indicates  whether  constraints  are  satisfied 

NTCTR    Number  of  the  point  on  which  the  program  is  working 

N405     Switch  showing  that  a  feasible  point  on  the  direction 
vector  could  not  be  found  and  the  negative  was  tried 

PO       Current  penalty  function  value   P(X) 

PX1      Penalty  function  value  of  lower  step  length  P(X3) 
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PX2  Penalty  function  value  of  upper  step  length 

RJ  Vector  of  previous  values  of  the  constraints 

TAL  Current  step  length   0 

TALI  Lower  step  length  0T 

Li 

TAL2  Upper  step  length  0 

TINC  Increment  by  which  step  length  increased  or 
decreased 

X  Current  X  vector  XQ+  0s 

XX  Temporary  storage  used  for  switching  values 

XI  Initial  starting  point  vector  XQ 

X2  X  vector  at  upper  step  length  XQ+0  s 

X3  X  vector  of  lower  step  length   XQ+0Ls 
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ALONG  A  GIVEN  SEAQCH 


M  N  t  N  F 1 


RSIGMA,  RJ(20Q),  RFC 

,  DELXC(IOO),  RHOINfPJTIOt  EPSI, 

,  X2(100)t  X3(100),  X«2(i0C>» 

PGRAC(IOG),  CHGIlOOli 
NPHASE,  NSATIS 


SLEROUTINE  OPT 
C  ThIS  CPT  SUBROUTINE  FINDS  THE  MINIMUN 
C  VECTCF  USING  CUAORATIC  INTERPOLATION 

INFLICIT  REAL*8(A-H.0-Z) 

CCMMON/SHARE/  XQQO),  DELUOO),  A ( 100, 100  )  , N , V , 
IthMl 

CCNNCN  /VALUE/  F  ,  G , PQ  , 

CCMMON/CRST/  DELX(IOO) 
1TI-ETA0, 

2  PSIG1,  Gl,     XK100) 
3XP1(1C0),PR1, 

4  FF2, F 1 ,  Fl,  RJl(200)t  DOTT, 

5  FPEV3,A0ELX,  NTCTR,  NUMINIf 
AES(CUyMY)=DABS(DLMMY) 
N4C5=1 
CCTT=C. 
CC  2  K=1,N 

2      CCTT=CCTT+DELX(K)*CELXO(K) 
GC  TO  5 
2  CC  4  1=1, N 

4  CElXm=-DELX(I) 

5  CONTINUE 
C  INITIALIZATION 

Nf^NN  +  1 

NTCTR=NTCTR+1 

TAL=0.0 

TALl^O.O 

C*2. 

KTER=0 

NPEC=0 

TINC=1.0 

CC  7  K=1,N 

X1(K)=X(K) 

X3(K)=X(KI 
7      FX1=P0 

NFIRS=0 
C  STAFT  BRACKETING  PROCEDLRE 

7AL=TAL+TINC 
10  CC    13    K  =  1.N 

12  X<K)=X1(K)+DELX(K)*TAL 

K7ER=KTER+1 

CALL  EVALU 

GC  TO  (170,20,20),  NSATIS 
C  RECUCE  STEP  SIZE 
20     C=C5 

IF(NREC,GT.20)  GO  TO  280 

NPEC=NRED+1 

TINC=(TAL-TAL1)*D 

TAL=TAL-TINC 

GC  TO  1C 
C  HAVE  2  CR  MCRE  FEASIBLE  POINTS 


eEEN  FCLND 


30 

21 
22 


IF(NFIRS.EC.l) 

NFIRS-1 

CC  22  K  =  1,N 

X2(K)=X(K) 

FX2=PG 

TAL2=TAL 

IF(FX2.GE.PX1) 
C  INCREASE  STEP  SIZE 
25     TINC=TINC*D 

TAL=TAL+TINC 

NPEC=0 

GC  TO  10 
40  CONTINUE 
C  RECRCER  POINTS 


GO  TC  40 


GO    TO  20 


IF(TAL.LT 

IF(F0.L7 

XX=TAL 

TAL=TAL2 

TAL2=XX 

XX=FO 

FC=FX2 


TAL2)  GO  TO  54 
PX2)  GO  TO  50 
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45 
46 

47 

50 

51 


F>2=XX 

CC  45  K=1,N 

XX=X(K) 

X(K)=X2(K) 

X2(K)=XX 

GC  TO  54 

T*L1=TAL 

PX1=PC 

CC  47  K  =  1,N 

X3(K)»X(K) 

TH=TAL2 

GC  TO  35 

T/SL1  =  TAL2 

F>1=PX2 

CC  51  K  =  1,N 

X2(K)=X2(K) 

X2(K)*X(K) 

TAL2=TAL 

FX2=P0 
C  CHECK  STOPPING  CRITERION  OR  IF  SEARCH  DIRECTION  SH3ULC  eE 
C  REVERSEC 

IF(ABS(ABS(P0/PXl»-l.).LT.l.E-7)  GC  TO  220 

CC  52  K«1,N 

IF(ABS(ABS(X1(K)/X(K)  )-l.).GT.l.E-7)  GO  TC  52 
52      CCNTINLE 

GC  TO  2E0 
52     CCNTINLE 

GC  TO  3  5 
C  CHECK  STOPPING  CRITERION  OR  IF  SEARCH  DIRECTION  SHOULC  E5 
C  REVERSEC 

54  IF(AeS(ABS(P0/PXl)-l.).LT.l.E-7)  GO  TO  220 
CC  55  K=1,N 

IF(ABS(ABS(X1(K)/X(K)  )-l. J.GT.l.E-7)  GO  TC  56 

55  CONTINUE 
GC  TO  2SG 
CCNTINLE 

IF(KTER.GE.IOO)  GO  TO  280 
IF(F0.GT.PX2)  GO  TO  46 
IF(PO.GE.PXl)  GO  TO  21 
CC  57  K=1,N 

IF(ABS(ABS(X2(K)/X(K)  1-1. ) .G£. l.E-7 )  GO  TC 
CCNTINLE 
GC  TC  2  20 

C  CLADRA7IC  INTERPOLATION  BEGINS***************************' 

60  KTR=0 

61  K 1 P  =KTR  +1 
AL=G.5*((TAL**2-TAL2**2)*PX1+ (TA L 2**2-7 AL 1**2 ) *P0+ 

1(TAL1**2-TAL**2 

2)*FX2)/ ((TAL-TAL2 )*PX  1+ <T AL 2-TAL 1 )*P0+< T AL 1-TAL > *P X2 ) 

IF(TAL.GT.AL)  GO  TC  75 
C  THROW  AWAY  LEFT  SIDE 

TAL1=TAL 

FX1=P0 

CC  71  K=1,N 
71     X2(K)*X(K) 

GC  TO  78 
C  THPCh  AWAY  RIGhT  SICE 
75 


56 


57 


60 


77 

78 


80 


TAL2=TAL 

PX2=P0 

CC  77  K  =  1,N 

X2(K)=X(K) 

TAL=AL 

CC  79  K^l.N 

X(K)=X1(K)+DELX(K)*TAL 

CALL  EVALU 

GC    TO    (170f80t85)t    NSATTS 

DC  82  K=lfN 


C  CHECK  STOPPING  CRITERION  OR  IF  SEARCH  DIRECTION  SHOULD  E  = 

GO  TC  82 
82 


IF(AES(ABS(X2(K)/X(K) ) -1. ) . GE. 1 . E-7 ) 

CONTINUE 

GC  TO  220 
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83 
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IFUES(4BS(P0/PX1  )- 

IF(KTR.GT.20)    GO    TO 

GC    TO    6  1 

NREC=0 

KTER=0 

GC    TO    20 
17C    CC    180    K=l,y 

IF(RJ(K))     190,190,180 
160    CCNTINUE 

NSATIS=4 

RETURN 
190  MN=C 

CC  200  K=1,M 
2CC  PJ1(K)=RJ(K) 

RETURN 
22C  CCNTINUE 
RETURN  hITH  SMALLEST  VALUE 

IF(FXl.GE.PO)  GO  TC  240 

PC=FX1 

CC  230  K=1,N 

X(K)=X3  (K) 

CCNTINUE 

RETLRN 
ERSE  SEARCH  0IR5CTICNS 

GC  TO  (290,300),N405 

N4G5=2 

NTCTR=NTCTR-1 

GC  TC  3 
fcPITE(6,580) 
CALL  TINEC 
C4LL  GUTPUT(l) 
C*LL  REJECT 

FCFMATC  80H   OPT  CAN-T 
1GIVES  A  LOWER  VALUE  OF 
STCP  22C42 
ENC 


)-lo ).LT„l.E-7) 
rn  220 


GG  TO  220 


230 
240 

C  RE\ 
2EC 
29C 


2CC 


5EC 


FIND  A  .FEASIBLE 
THE  P-FUNCTIGN. 


PCINT,ThAT 
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